\(\left(\mathrm{RC} 1^{\mathrm{st}} \text { ed: } 2.5\right)\) The RANDU generator, once popular on IBM machines, is based on the recursion: \[ X_{n+1}=65539 X_{n} \quad \bmod 2^{31} \] Illustrate the undesirable behavior of this generator with a computer experiment, and produce a \(3 \mathrm{D}\) plot which demonstrates the problem clearly. (Hint: show \(X_{t+1}=\left(6 X_{t}-9 X_{t-1}\right)\) mod \(2^{31}\)
Prove:
\[ \begin{aligned} X_{n+1} &=65539 X_{n}=\left(2^{16}-3+6\right) X_{n} \\ \left(2^{16}-3\right) X_{n} &=\left(2^{16}-3\right)\left(2^{16}+3\right) X_{n-1} \\ &=\left(2^{32}-9\right) X_{n-1} \\ &=-9 X_{n-1} \bmod 2^{31} \\ \Rightarrow X_{n+1} &=\left(6 X_{n}-9 X_{n-1}\right) \bmod 2^{31} \end{aligned} \]
\(2^{31} = 65339 | X_{n+1}-6 X_{n}+9 X_{n-1}\)
Now we may proceed into making a plot:
x = rep(0,1000)
y = rep(0,1000)
z = rep(0,1000)
for(i in 1:1000){
x[i] = sample(1:2^(31), 1)/(2^31)
y[i] = (65539 * x[i]) %% 1
z[i] = (6*y[i] - 9*x[i]) %% (1)
}
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.0001)
In class we discussed the Box-Muller algorithm for generating pairs of iid \(N(0,1)\) random variates.
a Show that if \(X_{1}\) and \(X_{2}\) are independent standard normal random variables then their polar coordinates \(R=\sqrt{X_{1}^{2}+X_{2}^{2}}\) and \(\theta=\) \(\tan ^{-1} X_{1} / X_{2}\) are also independent and derive their distributions.
b Show how to generate \(R\) and \(\theta\) by the inverse cdf method, and a pair of standard normals via the substitution \(X_{1}=R \cos \theta\) and \(X_{2}=\) \(R \sin \theta\)
c Implement the Box-Muller algorithm, and compare with the method described in class based on \(k=12\) independent uniform \([0,1]\) r.v. \(\mathrm{s}\) and the CLT. Are the two “statistically indistinguishable”? How do the respective computational costs compare?
Prove:
So \[ \begin{aligned} R^{2} & \sim \operatorname{Exp}\left(\frac{1}{2}\right) \\ \theta & \sim U[0,2 \pi] \end{aligned} \]
Sample \(U_{1} \sim U_{1}[0,1]\) and let \(R=\sqrt{-2 \log U_{1}}\), sample \(U_{2} \sim U[0,1]\) and let \(\theta=2 \pi U\) and \[ \begin{array}{l} {X_{1}=\sqrt{-2 \log U_{1}} \cos \left(2 \pi U_{2}\right)} \\ {X_{2}=\sqrt{-2 \log U_{2}} \sin \left(2 \pi U_{2}\right)} \end{array} \]
n = 1000
bm = rep(0,n)
for(i in 1:(n/2)){
U1 = runif(1)
U2 = runif(1)
R = sqrt(-2*log(U1))
theta = 2*pi*U2
bm[2*i-1] = R*cos(theta)
bm[2*i] = R*sin(theta)
}
## clt
clt = rep(0,n)
for(i in 1:n){
U = runif(12)
clt[i] = sum(U)-6
}
ks.test(bm, clt)
##
## Two-sample Kolmogorov-Smirnov test
##
## data: bm and clt
## D = 0.042, p-value = 0.341
## alternative hypothesis: two-sided
ks.test(bm, pnorm)
##
## One-sample Kolmogorov-Smirnov test
##
## data: bm
## D = 0.031265, p-value = 0.2824
## alternative hypothesis: two-sided
ks.test(clt, pnorm)
##
## One-sample Kolmogorov-Smirnov test
##
## data: clt
## D = 0.023099, p-value = 0.6601
## alternative hypothesis: two-sided
To generate the Beta distribution \(\mathcal{B} e(\alpha, \beta)\) we can use the following representation:
Proof:
random_beta = function(n, alpha, beta){
X = rep(0, n)
for(i in 1:n){
Y1 = rgamma(1, alpha, 1)
Y2 = rgamma(1, beta, 1)
X[i] = Y1/(Y1+Y2)
}
return(X)
}
reject_sampling_beta = function(n, alpha, beta, C){
X = rep(0, n)
j = 1
for(i in 1:n){
Y = runif(1)
U = runif(1)
if(U < dbeta(Y, alpha, beta)/(C*1)){
X[j] = Y
j = j + 1
}
}
return (X[1:j])
}
## alpha = 1.5, beta = 4
rb = random_beta(100, 1.5, 4)
rs = reject_sampling_beta(200, 1.5, 4, 2)
ks.test(rb, function(x){return(pbeta(x,1.5,4))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rb
## D = 0.093788, p-value = 0.3426
## alternative hypothesis: two-sided
ks.test(rs, function(x){return(pbeta(x,1.5,4))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rs
## D = 0.091995, p-value = 0.4069
## alternative hypothesis: two-sided
## alpha = 2, beta = 2
rb = random_beta(100, 2, 2)
rs = reject_sampling_beta(200, 2, 2, 2)
ks.test(rb, function(x){return(pbeta(x,2,2))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rb
## D = 0.10715, p-value = 0.2011
## alternative hypothesis: two-sided
ks.test(rs, function(x){return(pbeta(x,2,2))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rs
## D = 0.044987, p-value = 0.9782
## alternative hypothesis: two-sided
## alpha = 4, beta = 1.5
rb = random_beta(100, 4, 1.5)
rs = reject_sampling_beta(200, 4, 1.5, 2)
ks.test(rb, function(x){return(pbeta(x,4,1.5))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rb
## D = 0.12751, p-value = 0.07739
## alternative hypothesis: two-sided
ks.test(rs, function(x){return(pbeta(x,4,1.5))})
##
## One-sample Kolmogorov-Smirnov test
##
## data: rs
## D = 0.055482, p-value = 0.904
## alternative hypothesis: two-sided
This shows that the sampling method is still effective. But, depending on the value of c, the Accept-Rejection sampling algorithm is less efficient depending on the value.
For the Accept-Reject algorithm [A.4], with \(f\) and \(g\) properly normalized,
Proof
Without violating the assumption that \(f(x) < Mg(x)\), we may push M to its limit such that \[ M = \max \frac{f}{g} \] In this way, we end up with that the algorthm is much more effcient. If \(M\) is hard to calculate the exact form, we may sacrifice a bit by taking M is bigger than enough and ensure the correctness.
This is different from the given formula.
Given a normal distribution \(N(0,1)\) restricted to \(\mathbb{R}^{+},\) construct an AcceptReject algorithm based on \(E x p(\lambda)\) and optimize in \(\lambda .\) Repeat for intervals \([1, \infty),[2, \infty),[4, \infty) .\) For each interval, compare the efficiency to direct rejection sampling from the normal.
Prove:
\[\begin{aligned} M(\lambda) &=\max \left(\frac{f(x)}{g(x)}\right) \\ &=\max \left(\frac{\frac{1}{\sqrt{2 \pi}(1-\Phi(a))} e^{-\frac{x^{2}}{2}}}{\lambda e^{-\lambda x}}\right) \end{aligned}\] \[\begin{aligned} \frac{\partial \log M(\lambda)}{\lambda}=\left\{\begin{array}{ll} {\lambda-\frac{1}{\lambda}=0, \lambda \geq a} & {\Rightarrow \lambda=1} \\ {a-\frac{1}{\lambda}, \lambda \leq a} & {\Rightarrow \lambda=\frac{1}{a}} \end{array}\right.\\ \lambda=\left\{\begin{array}{ll} {1} & {x \in[0, \infty)} \\ {1} & {x \in[1, \infty)} \\ {\frac{1}{2}} & {x \in[2, \infty)} \\ {\frac{1}{4}} & {x \in[4, \infty)} \end{array}\right. \end{aligned}\]In this way, we end up with
\[M=\left\{\begin{array}{ll}{\frac{2}{\sqrt{2 \pi}} e^{\frac{1}{2}}} & {x \in[0, \infty)} \\ {\frac{1}{\sqrt{2 \pi}(1-\Phi(1))} e^{\frac{1}{2}}} & {x \in[1, \infty)} \\ {\frac{2}{\sqrt{2 \pi}(1-\Phi(2))} e^{-1}} & {x \in[2, \infty)} \\ {\frac{4}{\sqrt{2 \pi}(1-\Phi(4))} e^{-7}} & {x \in[4, \infty)}\end{array}\right. \text{ and } P(\text { Accept })=\left\{\begin{array}{ll} {0.76} & {x \in[0, \infty)} \\ {0.24} & {x \in[1, \infty)} \\ {0.08} & {x \in[2, \infty)} \\ {0.02} & {x \in[4, \infty)} \end{array}\right.\]
From_Exp = function(n, lambda, start, M){
res = rep(0,n)
j = 0
for(i in 1:n){
Y = rexp(1, lambda)
U = runif(1)
if(Y > start && U < dnorm(Y)/((1-pnorm(start))*M*dexp(Y, lambda))){
j = j+1
res[j] = Y
}
}
return (res[1:j])
}
length(From_Exp(5000, 1, 0, 1/(sqrt(2*pi)*(1-pnorm(0)))*exp(0.5)))/5000
## [1] 0.7538
length(From_Exp(5000, 1, 1, 1/(sqrt(2*pi)*(1-pnorm(1)))*exp(0.5)))/5000
## [1] 0.2486
length(From_Exp(5000, 0.5, 2, 2/(sqrt(2*pi)*(1-pnorm(2)))*exp(-1)))/5000
## [1] 0.0758
length(From_Exp(5000, 0.25, 4, 4/(sqrt(2*pi)*(1-pnorm(4)))*exp(-7)))/5000
## [1] 0.0238
Reject = function(n, start){
res = rep(0,n)
j = 0
for(i in 1:n){
Y = rnorm(1)
if(Y>start){
j = j+1
res[j] = Y
}
}
if(j==0) return(NULL)
return (res[1:j])
}
## a = 0
length(Reject(5000, 0))/5000
## [1] 0.5016
length(Reject(5000, 1))/5000
## [1] 0.1594
length(Reject(5000, 2))/5000
## [1] 0.0264
length(Reject(5000, 4))/5000
## [1] 0